﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SmartMathLibrary;

namespace SmartMathLibrary
{
    /// <summary>
    /// This class provides an incomplete LU factorization for square matrices.
    /// </summary>
    [Serializable]
    public class IncompleteLUFactorization
    {
        /// <summary>
        /// This field saves the matrix for the incomplete LU factorization.
        /// </summary>
        private SquareMatrix sourceMatrix;

        /// <summary>
        /// Initializes a new instance of the <see cref="IncompleteLUFactorization"/> class.
        /// </summary>
        /// <param name="sourceMatrix">The matrix for the incomplete LU factorization.</param>
        public IncompleteLUFactorization(SquareMatrix sourceMatrix)
        {
            if (sourceMatrix == (SquareMatrix) null)
            {
                throw new ArgumentNullException("sourceMatrix");
            }

            this.sourceMatrix = sourceMatrix;
        }

        /// <summary>
        /// Gets or sets the matrix for the incomplete LU factorization.
        /// </summary>
        /// <value>The matrix for the incomplete LU factorization.</value>
        public SquareMatrix SourceMatrix
        {
            get { return sourceMatrix; }
            set { sourceMatrix = value; }
        }

        /// <summary>
        /// Executes the incomplete LU factorization of the specified matrix.
        /// </summary>
        /// <returns>The LU matrix of the specified matrix.</returns>
        public SquareMatrix ExecuteFactorization()
        {
            int n = this.sourceMatrix.Rows;

            SquareMatrix matrix = this.sourceMatrix;

            for (int r = 0; r < n; r++)
            {
                double d = 1 / matrix.GetValueAtPosition(r, r);

                for (int i = r + 1; i < n; i++)
                {
                    double e = d * matrix.GetValueAtPosition(i, r);

                    this.sourceMatrix.SetValueAtPosition(i, r, e);

                    for (int j = r + 1; j < n; j++)
                    {
                        matrix.SetValueAtPosition(i, j,
                                                  matrix.GetValueAtPosition(i, j) - e * matrix.GetValueAtPosition(r, j));
                    }
                }
            }

            return matrix;
        }
    }
}